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Abstract 

In a relativistic setting, hydrodynamic calculations which include shear viscosity (which is first 
order in an expansion in gradients of the flow velocity) are unstable and acausal unless they 
also include terms to second order in gradients. To date such terms have only been computed 
in supersymmetric A/'=4 Super- Yang-Mills theory at infinite coupling. Here we compute these 
second-order hydrodynamic coefficients in weakly coupled QCD, perturbatively to leading order in 
the QCD coupling, using kinetic theory. We also compute them in QED and scalar theory. 



I. INTRODUCTION AND RESULTS 



Recently the Relativistic Heavy Ion Collider (RHIC) at Brookhaven has successfully 
created the quark-gluon plasma. Measurements of elliptic flow [1] indicate collective fluid 
behavior which implies a startlingly low viscosity [2]. Actually, measured in Poise the 
viscosity is enormously large; but this is expected of such a hot and dense system. It has 
recently been argued [3] that viscosity naturally scales with entropy density. Their ratio r]/ s 
is dimensionless [in natural units, used throughout; restoring h and c, it has units of ti\ and 
is conjectured to be bounded below by r^/s > l/47r (see however [4]). 

It is believed that the quark-gluon plasma created at RHIC displays a viscosity relatively 
close to this bound. But it is important to quantify this by comparing experimental results 
for elliptic flow spectra to the predictions of viscous hydrodynamics simulations. Several 
groups are engaged in this [5-9], but it is not as simple as adding a viscosity term to the 
ideal hydrodynamical equations. Indeed, it has been known for decades that relativistic 
Navicr-Stokes equations are acausal and unstable [10-12].-'^ 

Viscosity is just the flrst-order term in a gradient-expansion of corrections to ideal Eu- 
lerian hydrodynamics; Israel and Stewart showed 30 years ago that the stability problems 
could be repaired by the inclusion of certain second-order terms as well [11]. This is the 
guiding philosophy for most recent viscous hydrodynamics studies of the quark-gluon plasma. 

However, once one allows for some second-order in gradients terms, it seems wise to at 
least consider all second-order terms which could appear and to make an estimate of their 
size relative to the shear viscosity. This program was begun recently by Baier et al [13], who 
showed that, with the additional assumption of conformal invariance (a good approximation 
in QCD if the temperature is well above the QCD transition/crossover temperature of ~ 170 
MeV), there are five second-order coefficients, one of which is only relevant in curved space. 

It would be valuable to have a reasonable estimate of the size of these second-order 
coefficients, or an estimate of how they scale with the shear viscosity. Baier et al and the 
Tata group [14] have given one estimate, by evaluating the five coefficients in a toy model for 
QCD, strongly coupled A/'=4 Super- Yang-Mills theory (see also [15]). Here we evaluate the 
five second-order coefficients in QCD to leading order in the weak coupling expansion, using 
kinetic theory. In the thermal field theory setting the coupling expansion is not believed 
to converge very well (sec for instance [16]), so weakly coupled QCD should also be viewed 
as a "toy model" for QCD at realistic couplings. However we hope that the combined 
insight from the two "toy models" give a reasonable idea of the expected scaling of these 
second-order coefficients relative to shear viscosity. 

We evaluate the flat-space coefficients in Section II and the curved-space coefficient in 
Section HI. We then give an extensive discussion, in Section IV, of the physical interpretation 
of each second-order transport coefficient, and some interesting physical issues which arise in 
their computation. Certain technical details involving nonlinear corrections arising through 
plasma screening are postponed to Appendix A. But we will finish introducing the problem 
and present the main results and conclusions here. 



The easy way to understand this is to note that Navicr-Stokes equations are Euler equations plus a 
momentum-diffusion term, with the viscosity as the momentum-diffusion coefficient. But diffusion equa- 
tions possess infinite propagation speeds for information, which is problematic in a relativistic setting. 
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All hydrodynamic approaches are based on stress-energy conservation, 



0, 



which is 4 equations for 10 unknowns.^ The other 6 equations are established by gradient 
expanding the form of T'^^ about its equilibrium form. In the absence of nonzero conserved 
charge densities (which we will assume henceforward), in equilibrium^ 



T"" = (e + P)uV + Pg"", = -1 with k° > , P = P(e) . 



[1.2) 



This determines T^^'^ in terms of 4 unknowns, the energy density e and 3 components of the 
flow 4-vector u'^. However, if e,u'^ vary in space and time^ then we expect corrections to 
Eq. (1.2). For slowly varying e and u^^ the corrections can be expanded in gradients of these 
quantities. At first order in gradients and in a conformal theory, defining the rest-frame 
spatial projector 

A^"" = g^"" + u^u" (1.3) 
and working in flat space (so = and g^u = ry^,/), the only possible combination is 

(1.4) 

Here r] = rj{e) is the shear viscosity, defined as the coefficient multiplying the traceless 
part of the transverse symmetrized shear flow tensor. The bulk viscosity, deflned as the 
proportionality constant for the pure-trace part U'^" oc A'^'^ Aapd^'u^ , vanishes in a conformal 
theory. 

Baier et al ([13] Eq. (3.11)) show that there are four possible second-order flat-space 
terms: 



^^2 order 



+A2 



+ Ai 



- -A^ aafjCT 



a/3 



3 

+ A3 



1 

— I 

3 



,(1.5) 



= iA^„A,^(9"^^ - df'u'^) [vorticity] 



Physically, rn tells how quickly the anisotropic stress U^^" relaxes to the leading-order form 
—Tja'^'^, if it starts out with a different value. The parameter Ai tells how nonlinear the 
viscous effects are; A2,3 are similar but for systems with nonzero vorticity. An additional 
term k,{R^'' -|- . . .) is possible in curved space. It is these quantities we want to determine in 
weakly coupled QCD. We describe their physical signiflcance in more detail in Section IV. 

Expressing 1] in terms of the dimensionless ratio rj/s disguises the fact that rj really reports 
a time scale, roughly speaking the equilibration time of the system. The gradient expansion 
of Eq. (1.4), Eq. (1.5) is an expansion in this time scale divided by the scale of spacetime 



^ Here we only consider systems with vanishing densities of other conserved charges such as baryon number. 

^ We use the [ — metric convention 

^ It is also necessary to choose some convention defining e and u. We take the Landau-Lifshitz convention 
that U|J,T^"' oc u"" and e = u^u^T^"" . 
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variation of the system. To identify the time scale, divide rj not by the entropy density but 
by the enthalpy density: time. In A/'=4 SYM theory at strong coupling the 

ratio is = In weakly coupled QCD it is parametrically ~ giT\n{i/g) [17-20]. 

Similarly, the ratio of each second order coefficient to (e+-P) yields the square of a time. It 
is natural to expect ^ ~ (i+p)^- "^^^ numerical value of the ratio j^p^/i-^y = ^^^^^^^ 
is a convenient way to express the relative size of the second-order coefficient Ai to 77. In 
particular we expect most coupling dependence to cancel in this ratio, which should therefore 
differ relatively little between weak and realistic coupling. 

We find that at weak coupling, at leading order the ratios of second-order to first-order 
hydrodynamic coefficients are 



{e+P)riTn 
rf 
{e+P)K 
rf 

{€+P)\l 



5.9 to 5.0 (varies with g) : 6.10517 in Xcf)"^ theory , (1.6) 
0, (1.7) 
5.2 to 4.1 (varies with g) : 6.13264 in theory , (1.8) 



^2 



0. (1.10) 



The detailed coupling dependence of the two independent nonzero coefficients, rjTji and Ai, 
are shown in Fig. 1 and Fig. 2. These figures display results for QCD with either or 3 fiavors 
of quarks, and for e+e^ QED at realistic coupling = 5.9664 and = 5.4156] 

as well as indicating the results for weakly coupled A0^ theory. We have expressed the 
results in terms of ui-d/T the ratio of Debye screening length and temperature, which proves 
convenient computationally and is the right quantity for parametrizing whether a coupling is 
strong or weak at finite temperature. Numerically, ag = (2/127r)(mD/T)^ in 3-fiavor QCD, 
= (l/47r)(mD/T)2 in 0-fiavor QCD, and a^^ = (3/47r) (mo/Tj^ in 1-fiavor QED. Further 
discussion on these results and their physical meaning is postponed to the discussion section. 
Section IV. 

We find an exact relation X2/f]Tu = —2, in agreement with [21]. This relation is an 
automatic consequence of ultra- relativistic (conformal) kinetic theory. However unlike [21] 
we do not find Ai = rjTu- This is because [21] fixes an Ansatz for the functional form of the 
departure from equilibrium and drops some contributions arising from the nonlinearity of 
the collision operator. We discuss this in more detail in what follows. However in practice 
Xi/f]Tu is relatively close to 1. We also find that /t = = A3 in QCD, in QED, in scalar 0^ 
theory, and indeed in any conformal theory described by kinetic theory. But this does not 
mean that these coefficients are strictly zero; it means that they first arise in the perturbative 
expansion at a higher order than T^rn and Ai do. That is, Ai oc T^/ {g^ ln^(l/5f)) -|- 0{T'^/g^); 
but K may only scale as, say, T'^/g'^ and it is therefore zero in a leading-order evaluation, 
which only finds the oc T'^/g^ coefficients. This is discussed more in Section HI. 

For comparison, combining the results of [13] and [14] , the same coefficients in jV=4 SYM 
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FIG. 1: Coupling dependence of the ratio (e+P)rn/r/. This ratio compares the relaxation time 
scale for n^,^, tjj, to the time scale implied by the viscosity rj. 
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FIG. 2: Coupling dependence of the nonlinearity parameter Ai, expressed as the dimensionless 
ratio (e+P)Ai/77^. As explained in Section IV, there is an unresolved uncertainty in these curves, 
but it is smaller than the line widths. 
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theory are 



{e+P)r]Tu 
rf 
{e+P)K 

rf 

{e+P)Xi 
rf 

{e+P)X2 
rf 

(6+P)A3 

rf 



4 - 21n(2) ~ 2.6137, (1.11) 

4, (1.12) 

2, (1.13) 

-4 ln(2) ~ -2.7726, (1.14) 

0. (1.15) 



After scaling by the viscosity as described, the second-order coefficient r/rn is about twice 
as large at weak coupling as at ultra-strong coupling. The relation between Tyrn and A2 
valid at weak coupling is violated at strong coupling, and the coefficient Ai is also about 
two times larger at weak than at strong coupling. It is reasonable to expect that, in QCD 
at realistic couplings, the dimensionless ratios will fall between the weak-coupled values and 
the (generally smaller) ultra-strong coupled SYM values. Certainly we expect the QCD 
values for these dimensionless ratios to be of the same order of magnitude as what we find 
in both theories, wherever QCD is relatively close to conformal (starting somewhat above 
Tc). However given the difference in detail between values in the two theories it is tough to 
be confident in the exact values for realistic QCD. 



II. KINETIC THEORY TO SECOND ORDER 
A. Kinetic theory setup 

We will not discuss the derivation of kinetic theory here; for a review see [22-26]. Kinetic 
theory can be used when each of several criteria apply: 

1. There are long-lived quasiparticles (spectral functions for relevant fields or composite 
operators have sharp quasiparticle peaks). 

2. The density matrix is adequately approximated by a Gaussian approximation, that 
is, by the two-point function. Further, the system varies slowly in space and time, 
so we may work in terms of a space and momentum dependent distribution function 
/"(x, p). Here a is a label which runs over all quasiparticle types (species, spin, color, 
particle/antiparticle) . (Note that x and p don't commute, but if the spatial variation 
is slow enough then we can neglect the commutator and treat them as continuous, 
independent variables.) 

3. The quasiparticles dominate the measurables of interest and the dynamics. 

All of these criteria hold for weakly coupled relativistic field theories, even gauge theories 
[25], if we are interested in the transport coefficients which appear in the hydrodynamical 
description just discussed. The validity of the kinetic approach has been verified (at leading 
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order) by explicit diagrammatic analysis both in scalar field theory [27] and in gauge theory 
[28-30]. 

The kinetic theory description describes the time evolution of the distribution function 
/"(x, p). This is determined by the Boltzmann equation. In covariant notation, it is^ 

2P''dJ''{x,p) = -C[f] "Collision operator" (2.1) 

[np)llro^)[i±f^m)] - [i±r(p)] i[[i±r{k,)]f^m)) . 



Here we have defined p° in terms of the on-shell condition — Ep = A/p^~+m^ = p (in 
a conformal theory m = up to 0{g'^) medium corrections, which we will neglect since we 
seek a leading-order treatment), and we have introduced the shorthand 

The lefthand side of Eq. (2.1) describes the free propagation of particles; the time rate of 
change of the occupancy Edtf is determined by the particles' motion times the spatial 
variation of the distribution function dif{x,p). The righthand side describes the change in 
occupancy due to collisions, which are approximated as spacetime-local (so all / on the RHS 
are evaluated at the point x). The first product of population functions represents the rate 
at which particles of momentum p are scattered out of that momentum state; [l±/(k')] 
is a Bose stimulation (+) or Pauh blocking (— ) final state factor. The second product of 
population functions is the rate for the reverse process, producing a particle of momentum 
p. In equilibrium and in the local rest frame, [l±fcq{k)] — f{k)e^/'^ and so the two terms 
cancel by energy conservation, ensuring detailed balance. 

The Boltzmann equation rests on several approximations, such as the separation of scales 
between the distance between collisions {0{l/g'^T) in gauge theories) and the physical size of 
collisions {0{l/gT)) or deBroglie wavelengths of excitations {0{1/T)). It is not clear how to 
incorporate systematic corrections to these approximations. It is also problematic to evaluate 
the collision operator to high order in the coupling; for instance in QCD we anticipate that 
nonperturbative magnetic physics causes scatterings suppressed only by relative to the 
dominant 2^2 scattering processes. Indeed, we will shortly encounter (weak) logarithmic 
dependence on this scale in the second-order calculation performed here. Therefore it is 
not clear whether or how the kinetic treatment can compute transport coefficients beyond 



^ We use capital letters P for 4-vectors, boldface p for 3-vector components, and p for the magnitude 
|p| of the 3-vector. The collision operator here differs by a factor of 2p^ from that in [18-20]. This 

normalization difference will disappear when we integrate J^, since this integral carries a factor l/2po 
absent in [19, 20]. The overall minus sign on C is chosen so that its linearized form acts on the departure 
from equilibrium 5 f as a positive definite operator. To see the full covariance of the Boltzmann eqiiation, 
think of /"(x,p) as a function of 4-momentum P but with support only on the forward light cone, 

rix,p) = 5{p^)6{p'')rix,p). 
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leading order^ in g"^. So we will not try. This excuses us to simplify the collision operator to 
include only 2^2 and effective f ^ 2 scattering processes; in QCD the relevant collision 
terms are presented in [32] . It also means that we can neglect the scale dependence of the 
QCD coupling (the (3 function). Therefore QCD behaves as a conformal theory/ and the 
analysis of Baier et al [13] is relevant. 

B. Order by order expansion 

Our goal is to solve the Boltzmann equation for the case of a near-equilibrium system with 
slowly varying energy and momentum density (e, P*), or equivalently their dual variables, 
the temperature T and flow velocity it*. We write f{x, p) as a formal series 

/(^,p) = /o + A/i + AV2 + ... (2.3) 

with A a parameter keeping track of the order in derivatives. The lefthand side of the 
Boltzmann equation, Eq. (2.1) has an explicit derivative so it starts at 0{\). Therefore /o 
is fixed by the condition C[/o] = 0. The solution is (note that u^P^ < 0; /9 = 1/T as usual) 

/o(x,p) = (exp(-/?K'^P^)Tl)-\ (5^(5{x), u^^u^{x),p'^p, (2.4) 

with ^ = — for bosons and + for fermions. At first order we have 

2P^a^/o = -Ci[/i], (2.5) 

where we use the notation Ci to mean that C[f \ is expanded to first order in /i, see Eq. (2.36). 
At the second order we will have 

2P'^9^/i = -Cn[/i] -Ci;A,J/i] -Ci[/2] , (2.6) 

where Cu is the collision operator expanded to quadratic order in /i, Ci[/2] is the collision 
operator expanded to first order in /2, and Ci-^Mi the collision operator expanded to first 

order in /i and with the scattering matrix element also expanded to first order in /i. In 
principle there could also be a term Ci.„2 [/i] accounting for the /i dependence of particle 
dispersion relations, but this will be higher order in the gauge coupling so we can ignore it 
in this leading-order perturbative treatment.^ 

It is not our goal to determine the second-order departure from equilibrium /2. Rather, 
we only need to determine its contribution to the stress-energy tensor, which at leading 
order in coupling is determined in terms of / by 

T^u{x) = Y, I '^Pi^PJi^^ P) ■ (2-7) 



^ Note that the first corrections to the calculations we present here actually arise at order g, not g^. However 
we believe that the 0{g) corrections can be computed within kinetic theory; indeed this has been done in 
a few cases [16, 31]. 

^ For simplicity we will consider only masslcss QCD. 

* Dispersion corrections are 0{g'^) effects for the p ^ T particles which dominate transport coefficients. 
They are additionally suppressed because /i is chosen to have vanishing yoo(p) moment, and at order g'^ 
only this moment contributes to dispersion corrections for hard particles. 
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In particular this will mean that we only need spherical harmonic number £ = 2 components 
of /2. However since /i appears repeatedly in the expression Eq. (2.6) determining /2, we 
need its detailed form. Therefore the first step is to solve the first order Boltzmann equation, 
which was done already in [20]. So we begin by summarizing those results in the current 
notation. 



C. First order solution 

Explicitly evaluating the lefthand side of Eq. (2.5), 

2P''dM-f3P ■ u) = -2/^(-/3P ■u)(^P-u P'^d^P + (iPi'P^d^u^) . (2.8) 

Note that /o is a decreasing function so /q is negative. It is convenient to work noncovariantly 
at some point x and in the instantaneous rest frame at that point, so = 0, m° = 1 (using 
Roman letters for spatial indices, for which we will not distinguish between covariant and 
contravariant). At the point x the LHS of Eq. (2.5) becomes 

2P^^dM-PP ■ = 2/o(/5^) [E'^dtP + Pi{Ed,l3 - PEdtUi) - PiP^f^d^u,^ . (2.9) 

Separating the spherical harmonic number i — 2 and i — (traceless and pure-trace) parts 
of the last term, 

2piPjdiUj = (^PiPj - (^diUj + djUi - '^dijdkUk^ + '^E^dkUk , (2.10) 

the i — contributions in Eq. (2.9) are 

2f^E\dtl3 - pdiUi/^) (2.11) 

while the It — 1 term is 

2f,Ep,{d,(3 - (3dtu;) . (2.12) 

Note that, away from equilibrium, the definitions of (3 and m* are not unique; they are 
related to our choice of how to separate /o and /i, which is also not unique. The most 
sensible convention (Landau-Lifshitz) is to require in the local rest frame (the frame where 
7"°* = 2 JpPW"(p) — 0) that the departure /i + /2 + ■ ■ ■ carry no energy or momentum, 
^alpP^P^fiip) ~ ^- That means choosing the (undetermined) time derivatives dt/3 and 
dfUi such that the moments of the i — 0,1 terms vanish. At first order, this requires 

dtP = ^diUi and dtUi = ^diP (2.13) 

in the instantaneous rest frame; in covariant language 

ui'd^P = (^A'^'^d^u, and A'^'^u^d^u^ = ^A'^'d^P . (2.14) 

This fixes the definitions of f3 and u at first order in A. We will need these first-order 
relationships in evaluating the second-order departure in what follows. It also turns out to 
ensure that Eq. (2.11) and Eq. (2.12) cancel identically. 
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This leaves the i — 2 (traceless tensor) component as the sole source for the first-order 
departure from equilibrium, 

'^^am (m- - ^) f - cm , (2.15) 

where (jj^ was introduced in Eq. (1.4). It does not really matter whether uij multiplies p^Pj 
or ■pi'Pj — SijE"^ /"i in Eq. (2.15) since Uij projects out the trace piece; the latter shows the 
correct angular behavior, the former is simpler to use in some cases. 

A detailed treatment of the operator Ci[/i] is given in [18-20, 27]. What is relevant here 
is that is a rotationally invariant, linear operator on /i considered as a function of 

3-momcntum p. Therefore the angular structure of /i must match that of the lefthand side; 
/i must be of form 

/i(p) = '^{p^P,-5.,E'm'x{p)^^-fMv), 

= ^P^P^P^xi-pu- P) (covariantly), (2.16) 

with x(p) a dimensionless function of /3 and p = —u^P'^ which remains to be determined. 
By factoring out the powers of /? so % is dimensionless we have ensured that it is a function 
only of the dimensionless product (3p and not (3 and p separately. The relation between our 
notation and that of [19, 20] (AMY) is x = "i2(— /o)Xamy- The departure from equilibrium 
X is generically proportional to — /q = /o[l±/o] and it is also convenient to define a version 
where this is has been factored out, x — x/{~fo) ^-iid Xij — Xij/{~fo)- Note that x and x 
will both be negative definite. 

It is convenient to factor out 0-4^/2 from both sides of Eq. (2.15) and to consider it as an 
equation on the vector space oi i = 2 tensor functions of 3-momentum p. Using the inner 
product 

{A\B)^ I A{p)B{p) (2.17) 
Jp 

we can define Sij — 2{piPj — SijE'^/3), in which case the first-order Boltzmann equation is 

Pfo\Sij)^e^\xij). (2.18) 

At least formally we can then write 

\Xij)-PC^'fo\Sij). (2.19) 

The procedure for performing this inversion is described in [19, 20] and here we will simply 

assume that this part of the problem is already solved. Note in particular that besides 
explicitly scaling as g"^, the operator Ci also depends logarithmically on the coupling g due 
to screening efi'ects; therefore in gauge theories x(p) is a nontrivial function of g, as is 
anything which functionally depends on 

The first-order correction to the stress tensor is^ 

njj,l order = {Sij \ Xlm) = —^{Sij \ (3Ci Vo I ^Im) ■ (2.20) 



^ Our —{2E)Ci ^/g equals C^^^ of [19, 20]; our measure is 1/2E and our Sij is 2E times the normalization 
used there. These powers of 2E cancel to make the treatments equivalent. 
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In evaluating this quantity the relation for integrating over global angles holding relative 
angles fixed, 

with P2{x) the second Legendre polynomial, is useful. 
D. Second order treatment 

Now we roll up our sleeves and continue to the next order. Returning to Eq. (2.6), we 
will find that, formally, 

h = -Cr' (2P'^a^/i + Cn + C^m) ■ (2-22) 

Therefore wc need to compute the three terms on the righthand side, treating the first- 
order departure from equilibrium xif^E) as already determined. Actually we only need to 
calculate that part of /2 which contributes to the off-diagonal stress tensor 

n2order = i^ij I 

= -{S,j\C^'\2P^^dJ,+Cn[fi]+Ci.MAfi])- (2-23) 

But 

{S,j I Cr' = {Xij I Wnr' = -T{Xij\ (2.24) 

is known; therefore we need 

nLder = T {xrj I '^P'd.h + Cn[/i] + C,,MAh]) ■ (2-25) 

In other words we need the p integral, weighted with Xiji of three terms. No new operator 
inversions are required, though evaluating Cn and Ci;Xi require performing complicated 
integrals. 

1. ^P^'d^h term 

We begin with the 2Pi^d^fi term. This contributes to the most coefficients (r^Tn, Ai, 
and A2) but is the most similar to what we have already encountered. Wc compute it by 
evaluating 2P^(9^/i directly, taking the integral moment only at the end (but feeling free to 
drop terms which will vanish on angular integration). 

Since we are taking its spacetime derivatives, it is necessary to use the covariant form for 
/i, Eq. (2.16). The derivative can act on cr^^, on /3, or on x's argument; 

P"5a(/9V^''^'x(-/9«-^)) = /3V^''^"X(-) x3P"a,,ln/3 + /?3p-P^P''(5„a^,)x(..) 

-P^a^,P''P''x'{-)P'"P^{u^da ln/3 + . (2.26) 

We only need terms which in the rest frame are even in p. In the first and third terms 
a^iy^s indices are spatial so P^ and P'^ must also be; therefore in the first term there is only 
a contribution from Edt(3 and in the third term there is a contribution —E^dtP and p^jPdiUj 
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since dtUt = at rest. The middle term is trickier; do(Jij can be nonzero but so can diaoj; aoj 
vanishes only at x = but varies from away from the origin (the rest frame at neighboring 
points is not the same as at the origin). We can re-express it using diaoj = u^diapj, and 

((7"%) = a^(0) = ^ upd^a^P ^-a^^d^up. (2.27) 

In other words, 

diCTQj = -GkjdiUk . (2.28) 
Therefore the terms even in spatial indices are (also using Eq. (2.13)) 

(2.29) 

In the second term, the quantity in parenthesis is p^p"^(Ji„j2. In the first term we need to 
rewrite djUk-, decomposing it into its traceless symmetric, antisymmetric, and trace compo- 
nents; 

djUk = + 

_ djUk + dkUj - 25jkdiUi/3 Sjkdiui djUk - dkUj 
2 3 2 

= ^ + nj, + ^SjkdiUi. (2.30) 
Therefore this first term turns into 

l3^xi--)p'p'E(dtaij + ^aijdkUk - aikOjk - luik^j^ . (2.31) 
This term's contribution to Iljj 2 order is 

nij,2 order D (dtaim + -ClmdkUk — CrikCrmk — ^aik^mk^ X 

jp {piP, - {piPm - mm . (2.32) 

Using Eq. (2.21) the angular integration gives 2p^/15, replacing the Im indices with ij, re- 
moving trace parts, and leaving the radial integral /5^(307r^)~^ J pdpp^x{p)x{p) the overall 
coefficient. This contributes (with negative coefficient) to Ai and is the sole contributor to 
the terms Tn and A2, fixing the relation A2 = —2r]Txi, regardless of the form of the collision 
operator (in agreement with Baier et al [13]). This relation seems to be a robust prediction 
of kinetic theory. ^'^ 

Similarly, the second term in Eq. (2.29) contributes (note that xx' < 0) 

nij,2 order D -2/3^ / XX ( PiPj " ( PlPm " ^-^^ ( p^p, - ^Al) ^'"'^'"^ . (2.33) 



This relation between A2 and Tn was long known [11] but always in the context of Grad's 14 moment 
method [33] ; we see here that it is independent of this particular approximation but is more general to 
ultrarelativistic kinetic theory. 



12 



Evaluating this requires a special case of the angular integration relation Eq. (2.44), which 
applied to this case gives 

f f P^^ij\f P^Slm\f p'^Srs\ CTimCTrs / Sijaira(Jlra\ 

A™ l"* " " —A"' - — j — = ['"■'''■ — —) ■ 

(2.34) 

This term contributes positively to Ai, and is about twice as large as the negative contri- 
bution from the first term; indeed if x is constant, then this factor of 2 is exact. Previous 
work [18] often used the Ansatz that x is constant and it is not too far from the case. In 
general, if the detailed form of x is known then evaluating these terms is straightforward. 



2. Cuih] term 



Now consider The specific form of the collision operator now becomes relevant; 

we will first consider the case of a 2 <-> 2 coUision operator. It is convenient [18] to introduce 

M(3E) = -f^{(3E)M(3E) (2.35) 

and similarly for /2. Writing / = /o — /o(/i + h), we find to second order, 

/(p)/(k)[l±/(p')][l±/(k')] - [l±/(p)][l±/(k)]/(p')/(k') 
= Mp)fo{k)[l±fo{p'm±fo{k')] X 

([A(p)+/i(k)-/i(p')-/i(k')] + [/2(p)+/2(k)-/2(p')-/2(k')" 

+ /i(p)/i(k)/o(p)/o(fc)(e^ - 1) + A(p')/i(k')/o(p')/o(A:')(l - e"^) 

+ [/i(p)/i(p')/o(p)/o(p') (e* - e^) + {p' ^ k') + {p^k) + {p,p' ^ k, k')] )(2.36) 

plus terms which are third order in gradients. Here (p' — > k') means the first term in the 
square brackets, but with the substitution p' —>■ k'. The first two square-bracketed terms 
are responsible for Ci[/i] and Ci[/2]; the last two lines are quadratic in fi and are therefore 
what we meant by Cn terms. The contribution of Cn to Uij will involve 

n,,- 2 order D / {27r)'6\P+K - P' - K')\M\' Mp) m)[l±Mp'm±Mk')] X (2.37) 

J pkp'k' 

TXij{p)^^^^Y^ Xim{p)Xrs{^)foip)fo{k) (e^ -1^+5 more terms] . 
For collincar effective 1^2 processes wc similarly need {p' + k' = p) 

/(p)[l±/(p')][l±/(k')] - [l±/(p)]/(p')/(k') 
= /o(p)[l±/o(p')][l±/o(A:')]x (2.38) 



([/i(p)-/i(p')-/i(k') 



+ 



/2(p)-/2(p')-/2(k') 

+/i(p')/i(k')/o(p')/o(^')(l - e^) + [/i(p)/i(pO/o(p)/o(p') (e^ - e^) + {p' ^ k') 

The contribution to Ilj^ is of similar form to Eq. (2.37). These terms clearly depend in 
detail on the available processes and their matrix elements |A1p; they also require multi- 
dimensional integration over the external particle momenta. However the relevant matrix 
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elements and useful parameterizations for the angular integrations have already appeared 
[20], so we will concentrate on what is new, which is the angular structure. 

In evaluating Eq. (2.37) we will encounter an integration over global angles, keeping 
relative angles between p, p', k, k' fixed. Since the matrix elements do not depend on global 
angles, we may perform this global angular integration first. Introducing the notation 

for the traceless symmetrized part, the generic integral we need is of form 

^hH^ / P(iPj)k(ik„)p'^^p;^ , (2.40) 

where we will normalize so that f.o 1 = 1. We show how to deal with a slight general- 
ization of this form, needed in evaluating Ci-Mi ■ Consider 
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/ AijBi^Crs , A,B,C of form Aij = p^^q,) (2.41) 

^global 



that is, each A, B,C is a distinct traceless symmetric tensor. The global angular integration 
over AijBimCrs must give a rank-6 tensor, symmetric and traceless on each pair of indices. 
There is only one such tensor: 

f AijBi^Crs = C[A, B, C] (Sii6jrSms + 7 permut. - ^{SrsSuSjm + 5 permut.) 

•''f^global V 3 

+ Y^^AmSrs^ ■ (2.42) 

The coefficient C[A, B, C] is determined by contracting each side with SuSjrSms, yielding 

C[A, B, C] = ^AijBi^Cjm . (2.43) 



Therefore 



'—^ I AijBimCrs — — (cTiiaji ■^CTimO'lrn] ArgBrtCgt ■ (2.44) 

1 Jn^iob.i 35 V 3 J 

In particular, in evaluating Eq. (2.37) we will need angular moments of form 



O'lm.O'r 
4 



P(iPj>k(;k„)P(^p',) 

^global 



1 f 

35 \ 3 



'^Im'^lm^ ^Xpk^pp'^kp' ^pk -^pp' -^kp' (2-45) 



where we define Xp^ = p • k. This result together with results in [20] are sufficient to 
compute the Cn contribution. Note that the contraction of a tensors above is precisely the 
one defining the coefficient Ai in Eq. (1.5). Therefore the term Cn strictly contributes to Ai. 
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3- Ci-Mi contribution 

If we calculated x in a gauge theory, using the vacuum matrix elements, we would find 
a log divergence in Ci due to the Coulomb singularity, and therefore x would be zero. 
Therefore it is essential in applying kinetic theory in a gauge setting to include the physics 
of dynamical screening [18], both for gauge boson and for fermion exchange. 

However, dynamical screening depends on the density of plasma particles and their 
momentum distribution; the matrix element Ai is itself a function of /, Since 
/ = /o + /i + . . ., we can expand the matrix element as well; 

M[f] = M[fo] + A^/iW^lly + 0{X') (2.46) 

where as before A keeps track of orders in gradients. As shown in Eq. (2.36), the product of 
population functions in the collision operator is only nonzero at 0{X); therefore the 0{X) 
correction to first gives rise to a nonzero effect at second order in A. In particular 

C,.,MAfi] = J^^^^pT^m. . .) Jj,{r) (-^[/ol^^^ + h.c.) (2.47) 

x/o(p)/o(^)[l±/o(p')][l±/o(^')](/i(p) + /i(k) -/i(p') -/i(k')) ■ 

The contribution to Ilj^ is J Txij of this. 

The functional form of dAi / df is somewhat complicated but is only significant for small 
exchange momenta, that is, when one of the Mandelstam variables is small, say, t ^ m^- 
Therefore, in the context of a perturbative treatment it is fair to work in the small exchange 
momentum approximation, \t\ <^ s. This simplifies both the form of dM./df and of the 
integration structure. However the specific details for evaluating Ci-^Mi complicated 
enough that we have postponed them to Appendix A. 



III. KUBO FORMULA FOR m AND k 

The previous discussion has determined all but one of the second-order hydrodynamic 
coefficients; since we worked in flat space we were unable to determine the coefficient k. 
Here we evaluate k without leaving flat space, and provide an alternative evaluation of rn, 
by making use of a Kubo relation derived by Baier et al [13]. There it is shown that the 
two "linear" second-order coefficients, ryrn and k, can be determined if one can evaluate the 
retarded Green function for the stress tensor^^ 

^^-^-(a;, k) = y" d4xe-^'^*+^'^-^e(i) Tr pr [7;,(0) , 7;,(x)] (3.1) 

(with Pt the equilibrium, thermal density matrix) and expand it to second order in uj ^ at 
vanishing kx,ky. In particular (Eq. (3.14) of [13] in our conventions) 

Q'^yT.y(^^^ ^^-j ^ _-p + + ^ {u'^irim - k/2) - klK/2) . (3.2) 



Our convention for the retarded function is missing a factor of i found in many definitions; our retarded 
function for a free particle is Gr(P) = -i/{P'^ + + zep°) or Gr(P) = i/{p° - E + ie){p° + E + ie). 
Therefore 2Disc Gii{u)) = p{u)) the spectral function is real. 
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Note that all correlation functions in this section are for a plasma in equilibrium. 

We can use kinetic theory to compute a related equilibrium correlator, the Wightman 
function 

G>,T.yT.y(^i^^ k) = j d'^xe-^*+'^-'' Tr pTT^y{0)T^y{x) . (3.3) 
The relation between these correlation functions is that 

G>{uj, k) = ^—-{Gr{uj + le) - Gr{uj - ^e)) ~ -2 Re Gr{uj + te) . (3.4) 

(In the second relation we made the approximation uj T, valid for all frequencies of 
relevance here.) This relation can be inverted into a Kramers- Kronig relation 

To evaluate the Wightman function G^ , recall that the Fermi/Bose distributions have 
fluctuations which are independent for each a, p and of magnitude Sf{p) — /o[l±/o] = 
—fQ{p). The instantaneous value of T^y is 

Txy{x, t)^2 PxPySfip, X, t) , (3.6) 

Jp 

which averages to zero. But the two-point function does not; 

G>{x,t) = {T,y{0,0)T,y{x,t)) =4 [ p.PyP',p'y{6f{p',0,0)6f{p,x,t)) . (3.7) 



We can evaluate this at positive t by pretending that p^p^/o[l±/o](p') is a source for depar- 
ture from equilibrium in the Boltzmann equation and evaluating the expectation value for 
Txy with the resulting departure linearized^^ Sf{p,x,t). The relevant Boltzmann equation 
is 

^PxPyfMm^'ix) + 2{Edt+pA)Sh{p,x,t) = -Ci[6f{x,t)] . (3.8) 

The spatial Fourier transform is trivial, removing S^{x) and replacing di with iki. The time 
transform is more subtle. If C were replaced by a relaxation time C[f] — > 2EVfi and ignoring 
kz for the moment, we would have 

(5/[relax-time-approx;t] = -^^{-fo)2pxPy, 

(5/[relax-time-approx;a;] = (^ 2Er + 2iujE ^ ^"^") '^"'^"^ ^^'"^^ " ^^'^^ 

Instead C is an operator. Moving the spacetime derivatives to the righthand side and 
formally inverting, one finds 



(3.10) 



^"^ Sf is not quite the same as /i in the previous section; it includes terms second order in gradients but first 
order in the departure from equihbrium, that is, it wiU contain terms quadratic and higher in spacetime 
derivatives but is linear in m. 
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The stress-stress correlator is the value of T^y arising from this /i, which is 

G>M = (S,|V> = (S,| + -c.) /;|S.,> . (3.11) 

Now we use and the Kramers-Kronig relation to evaluate Gr. First consider the case 
where kz = but we allow u' to be finite. Then (combining fractions) 

Because C has a purely real and positive spectrum, we are free to perform the cu integral by 
the method of residues, enclosing only the pole arising from (C — i2Euj)] 

n=0 

(3.13) 

The leading term in the expansion is 

-z(S',,|^(-/o)|<S',,) ^-^2^1 (27r)3T p2 /o[l^/o] = "^^^s"^ ^^-l^) 

a 

which is I of the expected — iP. (Here — X]a(-'^(^-"^^'^^) '-"^ 1 (fcrmion)).) The remaining 
I of —iP arises from ~ T (large frequency cut) contributions to which we have not 
computed here, and which give only order contributions to r], rjTu, n, which we therefore 
neglect. 

The first subleading cx u' term in Eq. (3.13) reproduces Eq. (2.20) and the last term 
allows us to calculate the combination {rjTu ~ i^f^)'- 

Vru-K^ P{S,y I C^\2E)C^\-f;,) \ S,y) = T{x,y I 2E{-Q \ x.y) (3.15) 

which leads to the same result we had for jyrn previously in Eq. (2.32). This already shows 
us that K = 0. 

To establish that k = in another way, we directly evaluate the second order in k term 
at vanishing uj' . The retarded Green function is 

G^{u' = Q^kz) = (^.,1 / ^ ^ (- , + c.c.) {-f,)\S^y) . 

J 2t:T uj — uj' — le \Ci — i2Euj + ^2p^fc^ / 

(3.16) 

The ratio ujjiuj — uj' — ie) cancels. Because C has positive definite spectrum we can again 
perform the u integral by closing the contour above for the 1/(C — iEuj + ipzkz) term and 
below for the 1 / {C + lEuj — ip zkz) term. There are no poles to pick up, but there is a nonzero 
contribution from the contour-closing arc because the integrand only falls as 1/uj. However 
this arises in the extreme large uj region where the finite operators C,Pz are subdominant 



The integrand needs to be regular at w = for this cancellation to work, otherwise the ie prescription is 
nontrivial. However the good properties of C\ ensure this is the case. 
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and can be dropped. Therefore we find a kz independent result. Equivalently, we could 
Taylor expand about small kzi 

C + i2Eu - 2ipzkz ~ C + i2Eu ^ ^ C + 2iEu (^^^^^^c + 2iEu) ^'^'^^^ 

and integrate term by term; on all but the first k^ independent term the integrand falls as 
l/cij^ or faster, and we may close the contour away from all poles and pick up no contribution. 

Therefore the expansion of G-r{ui, k) in powers of kz at vanishing ui shows no k dependence, 
and the second-order coefficient k vanishes. To clarify, the expansion in nonzero kz and uj 
will contain nonvanishing terms, of order ukl etc. It is only the kz dependent terms at a; = 
(or vanishing order in u) which vanish in kinetic theory. Note that we did not have to make 
any assumptions about the collision operator C to arrive at this conclusion, except that it is 
space-local and positive definite (the equilibrium ensemble is stable against perturbations). 

This result is not too surprising. As explained in [13], another way of interpreting the k^ 
coefficient is that it gives the correction to the stress tensor if there is a spatially varying but 
time-independent traceless metric disturbance hxy{z) ^ 0. But examining classical phase- 
space trajectories for this specific background shows that an initially equilibrium distribution 
freely propagates to remain in equilibrium (at linearized order and when the geometry is 
time independent) . Explicitly, in curved space the Boltzmann equation is [34] 

p^d,,f{x,p) - r\,p''p''d^,f{x,p, t) = -C,[5f] . (3.18) 

For the case gf^u — rjf^i, + h^i,-! hxy — hyx — ae^^^ with a time independent and all other 
components zero, the nonzero Christoffel symbols are 

r r^xz r "^dx^ hxy ■ (3.19) 

Since F is already linear in h we may evaluate dp\f using the flat-space form for /q, dpxfo — 
{fo)P^/P- The second term on the lefthand side of Eq. (3.18) is therefore 

-r\,P^p''d,.fo = -^-^{Qdx^hxy . (3.20) 

To evaluate the first term, we have to evaluate / to first order in h. The equilibrium 
form is /o = {cxp{i3g^i,u^P'^) ^ 1)~^, and since only m° is nonzero and goi^ is unchanged this 
is /o = l/{e~^^ -F !)• However p° is defined implicitly in terms of via g^i.P'^P'^ = 0. 
Therefore p° = a/p^ -l- 2hxyP^py = p-\- hxyP^p^/p plus terms quadratic in h. Evaluating the 
space derivative therefore gives 

P^dx.fo = ^^^{Qdx^hxy . (3.21) 

The two terms cancel, meaning that the system remains exactly in equilibrium to linearized 
order in h. 

Since this argument relies only on classical phase space propagation, the coefficient k 
will first arise when this classical phase-space picture becomes insufficient. The parametric 
behavior of Ai ~ jg^ arose as T^/Z^fp, involving two powers of the mean free path. Our 
phase space argument shows that k must involve one power of the scale where classical phase 
space treatments break down, which is the scale set by the inverse dcBroglie wavelength 
T. Therefore we expect that k ~ T^/(/infpT) ~ T"^ / g* (at most). Computing the first 
nonvanishing contributions to n at weak coupling is beyond the scope of kinetic theory and 
of this work. 
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IV. DISCUSSION 



We clarify and discuss in turn the meaning and origin of the second-order coefficients 
within kinetic theory. In particular, consider shear flow with Ozz = ~2c, a^x = o^yy = c 
with c positive. This is Bjorken contraction, with some radial expansion to preserve volume 
(or pure Bjorken contraction plus a conformal transformation). In this case we expect a 
particle distribution to become prolate along the z axis, leading to T^z > T^x.Tyy. This 
is what happens. The magnitude, integrated over p^dp, determines rj. The deviation from 
equilibrium depends on p and is described by p^xip)^ the relative departure from equilibrium 
/i//o[l±/o] as a function of p. In a relaxation time approximation, x oc in a momentum 
diffusion approximation x 1- 

The physical meaning of rn is, how far Tzz comes from this expected form if the rate 
of Bjorken contraction is changing with time. If Bjorken contraction is speeding up, the 
particle distribution should reflect the smaller value which used to be valid; hence Tzz should 
be smaller, meaning the proportionality constant Tzz = —Wzz + V^ndtCzz should be positive 
(since azz is negative). This is the sign we obtain. But how much smaller? This depends on 
how quickly the distribution relaxes back to equilibrium. The size of 7]/{e+P) also depends 
on how quickly the distribution relaxes to equilibrium, so we expect some relationship rn ~ 
ri/{e+P). But the proportionality constant depends on whether all particles equilibrate in 
the same way, or some particles take longer to equilibrate. If high momentum particles 
take longer to relax to equilibrium, then they can store information about the value of 
azz further into the past. As a result, if we make a relaxation time approximation, then 
X ocl/p gives Tji = 5r]/{e+P), whereas the momentum diffusion approximation x cc 1 gives 
Tn = 6r]/{e+P). Figure 1 shows that the value moves from close to 6, at weak coupling, to 
nearly 5 at stronger coupling. This occurs because at weak coupling collisions are dominated 
by soft scattering, which acts like momentum diffusion and gives quite close to x oc 1 (see 
[19]), while at larger couplings coUinear splittings become more important and try to enforce 
X ocl/p (see [20]). So this coupling behavior is expected. 



y y 




Bjorken Contraction (Vorticity) rotation 

FIG. 3: Illustration of the physical origin of r] and of A2. Under Bjorken contraction (left), the 
momentum distribution becomes prolate long the z axis. But under rotation with dzVy > (right), 
the prolate axis gets rotated to have a y component, so Tyz > 0. 



The value in scalar A^'* theory is slightly higher than 6. However, if we replace Bose statistics with 
Boltzniann statistics, it turns out that the Ansatz % oc 1 is exact, and Tn(e+P)/77 = 6 exactly at leading 
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The relation between rn and A2, and the sign of A2, also have fairly simple interpretations. 
First the sign. Physically A2 tells what happens to a system which is both Bjorken contract- 
ing (nonzero a^z) and rotating (nonzero vorticity, say flzy > 0). As illustrated in Figure 3, 
in this case the contraction makes the particle distribution become prolate; but the vorticity 
skews this distribution so it is not aligned with the Bjorken contraction axis. That should 
lead to a positive Ty^, which for < and Qzy > requires A2 < 0. The proportionality 
constant depends on how large the original zz asymmetry was, which depends on 77, and 
on how long the induced xy skewed distribution "lives," which is set by rn. Accounting for 
numerical factors turns out to give A2 = —2'r)Tu, as we find. 

Next consider Ai. For our example of Bjorken contraction, 

^zz = -r](Jzz + Ai((T2/o-^i - 5^^o-f„/3) = 7]{2c) + Ai(2c^) . 

Therefore a positive Ai means that for Bjorken contraction, the stress tensor deviates further 
than normal from equilibrium. On the other hand, reversing the sign of c to consider 
Bjorken expansion, the deviation from the equilibrium value of Yl^z is reduced. Therefore 
Ai tells whether equilibration is accelerated for Bjorken expansion (Ai positive) or Bjorken 
contraction (Ai negative). 

Our calculation shows that there are 3 contributions to Ai. First, if the particle distri- 
bution has already become prolate, then further Bjorken contraction generates a different 
amount of prolateness than it would from a spherically symmetric distribution. This is 
the part contributed by 2P'^9,,/i. The sign turns out to be positive and the magnitude 
dominates all contributions to Ai. 

The contribution to Ai from Cu reflects the change, in going from a thermal to a prolate 
momentum distribution, in the set of scattering targets a particle has. Whether this accel- 
erates equilibration or slows it down depends on typical scattering angles in a rather compli- 
cated way, indicated by the rather complicated angular integrations involved in Eq. (2.45) 
and Eq. (2.37). This leads to considerable angular cancellation. For instance, in A0^ theory, 
where the matrix element = A^ shows no preference for particular scattering angles, the 
contribution to Ai from Cu is -1-0.0372. If we replace Bose with Boltzmann statistics in A0^ 
theory, the cancellation on angular averages becomes exact and Cu gives no contribution 
to Ai. In QCD the contribution is also small, due to significant angular cancellation; for 
3-fiavor QCD the Cu contribution to Ai varies between —0.18 at weak to —0.45 at stronger 
couphng. The negative sign means that prolate distributions show accelerated equilibration. 

The contribution to Ai from C-i-Mi reflects changes in the efficiency of scattering and 
collinear splitting because of changes in plasma screening. This is interesting because it is 
where the precursors of plasma instabilities (see [36-38]) can enter the game. An anisotropic 
particle distribution weakens the stabilizing effect of plasma screening for certain particle 
directions p and exchange momenta q. In particular, in directions where /i(p) is positive, 
these particles have enhanced scattering via soft magnetic (G? ) gluon exchange with q ± p. 
One might guess that this leads to a large negative contribution to Ai. However we find 
that extensive angular cancellations occur which make the contribution arising from elastic 



order in A. 

Ai does NOT indicate the "anomalous viscosity" expected from plasma instabilities [35]. "Anomalous 
viscosity," for which falls below the linear term for all flow patterns, would be indicated by a large 
positive value for the third order term Hij oc UijaimO'im- 
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scatterings very small, and free of IR divergences, see the discussion at the end of Appendix 
A3. 

The same does not happen for collinear splitting. If the particle distribution becomes 
prolate, the approach to equilibration would be accelerated (Ai < 0) if the particles traveling 
in the prolate (z) direction show a higher rate of collinear splitting, since such splitting is an 
equilibrating process. The rate of collinear splitting depends on the efficiency of transverse 
momentum diffusion. But the proto-plasma instability caused by a prolate distribution is 
automatically the right one to enhance such transverse momentum diffusion for particles 
moving along the z axis.^^ Therefore the contribution of collinear splitting processes in 
Ci-Mi should contribute negatively to Ai and give the first hints of the effects of plasma 
instabilities. 

The fractional change in scattering efficiency due to /i grows at small momentum ex- 
change as This behavior is expected; for weakly anisotropic plasmas only the smallest 
g's show plasma instabilities, which appear in perturbation theory to give an infinite scat- 
tering rate. Since is treated as formally infinitesimal, there is no finite momentum q 
which becomes unstable, but the restoring effect of the plasma is changed more and more 
for softer and softer magnetic q. This leads to an IR log divergence in the total momentum 
transfer rate / cPq±q'^C{q±), see Appendix A 5. Therefore the change to the rate of collinear 
splitting is log divergent when computed at leading perturbative order. 

This means that our result for Ai actually includes a (negative in sign) logarithmically 
divergent contribution, at least using the perturbative calculational tools we employ here. 
The log is ln(mD/e), with e an artificially imposed minimum momentum transfer, imple- 
mented by modifying q^ ^ q^ + e in the denominator for transverse gauge boson exchange 
when computing this process. 

Physically, there really will be a limit on the infrared end of momentum transfer. In 
QCD we expect e ~ g'^T the magnetic screening scale. This is where the perturbative 
treatment of plasma corrections to gauge field propagation breaks down. Unfortunately we 
cannot compute the exact form of this cutoff (the constant under the log, \n{mY) / g'^T) + k) 
because this momentum region is strongly coupled. Similarly, we expect that in QED 
the perturbative treatment of screening also breaks down for q ~ e^T, where the physical 
distance of particle propagation involved is of order the large-angle scattering length and the 
electron propagators cease to behave like Eikonal propagators (as assumed in the hard-loop 
computation of self-energies). It might be possible to compute the constant under the log, 
ln(mD/e*^T) -|- k , but we have not done so. 

As a result, we have not actually been able to compute the complete finite-coupling value 
of Ai. Rather, we have guessed what the cutoff e on transverse momentum should be; we set 
e = g'^T/2 in QCD and e = e'^T/10 in QED. This leaves an uncertainty in our results, set by 
the coefficient on the ln(mD/e) term arising from Ci-Mi from collinear splitting processes. 
Fortunately, it turns out that this contribution is numerically tiny. If the constant under 
the log shifts by 1 (the correct cutoff is g'^T/bA rather than g'^T/2) then our result for Ai 
changes by less than 0.003 in 3-fiavor QCD and less than 0.0003 in pure-glue QCD or QED. 



It also slows down transverse momentum diffusion for particles in the "cqiiator" of the prolate distribution, 
slowing their approach to equilibrium. The angle averaged rate of splitting remains constant at this order. 
But Ai does not depend on this angle averaged rate; it is dominated by what happens along the axis of 
prolateness (or oblateness). Therefore we can get the right sign by paying attention only to what happens 
to particles along the z axis. 
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The extreme smallness of this effect arises as the product of several small things. First, 
collinear splittings are not that important in driving thermalization. Second, the splitting 
rate is reduced in some directions, and there is some angular averaging which reduces the 
total importance of the shift in the splitting rate. Third, the change to the splitting rate in 
any specific direction also turns out to be numerically small. This is another indication that 
in practice the physical importance of plasma instabilities turns out not to be very large. 

We end the discussion by commenting about the range of validity of our calculation. In 
Figures 1, 2 we have plotted our results out to mn/T — 3, which corresponds to quite a 
large couphng ctg = .48 in 3-flavor and ctg = .72 in pure-glue QCD. The calculation certainly 
cannot be believed at such couplings; probably it becomes inadequate beyond m-o/T = 1 
(see [16] for a next-to-leading order calculation of a similar transport coefficient). The scaled 
results for rn and Ai are weakly dependent on details of the theory, as shown by the almost 
identical results for theory and QCD at weak and relatively strong couphng. But they 
rely in an essential way on the validity of kinetic theory. There will be 0{as) corrections 
which cannot be incorporated in kinetic theory, which we generically expect to change the 
shape of the curves and which we do not know how to compute. Therefore the flatness of 
the curves in the figures can only be taken seriously at small m-d/T (we would guess below 
mD/r= 1.5). 
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APPENDIX A: MATRIX ELEMENTS AT NONZERO a^^ 

A particle of momentum P scattering from a particle of momentum K via gauge boson 
exchange with a soft exchange momentum Q, \Q'^\ <^ \P ■ K\ does so with a leading order 
matrix element (suppressing group factors) 

M = 2P''G^,2K'' , G-^ = Q^g^, - Q^Q, - B^^/] + (Gauge fix) , (Al) 

with G, n understood as the retarded propagator and self-energy. What is relevant here is 
that n^,^ explicitly depends on the medium through its distribution function /. Write it as 
n/ii/f/] = n^j^eq + <^n^,^[/i] plus terms of higher order. Then the squared matrix element 
becomes 

|>1|2 = MoMl + MoM\ + MlMi + 0{\^) , 
Mo = 2P'^G'^,2K^ 

Ml = 2P'^G'^„5^"''[/l]G';3.2X^ (A2) 

where to simplify notation G now means the equilibrium propagator. Since B is suppressed 
relative to G~^ unless ~ g^T"^-, we can freely treat as small in what follows, system- 
atically expanding whenever possible in p. A; ^ 5, 5°. Similar expressions are also needed for 
fermionic exchange processes and the fermionic self-energy. 
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Our goal in this appendix is to evaluate Eq. (2.47). Clearly as a first step we need to 
evaluate 611°'^ and its fermionic equivalent; then we need to use this to evahiate (A^o-^i + 
h.c.) and perform the momentum integrations. In addition, the collinear splitting rate is 
sensitive to 511 because it depends on the rate of soft momentum exchange; so we will have 
to revisit the rate of coUinear splittings as well. 



1. Bosonic self-energy 

With the sign convention established in Eq. (Al), for soft 4-momentum Q = (5°,q) the 
leading order (retarded, hard-loop) self-energy is [39] 



d^p dfjp) 

Qpk 



v^g 7, 

V • q — — ie 



(A3) 



where v = p/p and p — |p| as usual. The sum is over species, spin and particle/ antiparticle 
but not color. Setting f — fo and using 



R 

recovers the usual HTL self-energies: in strict Coulomb gauge, which we use henceforth, 

q^ — UJ^ ^ TT lr\\ 2 { ^^ 



Goo= , ^ (n^ =^-^^^^ noo(g)=m^ 2--ln^ , (A5) 

-r - noo((5) r \ Q ^-q J 

^ij = T2 72 — "FTTT^ = - liljjGT , nT(Q) = -m — + — In 



(Throughout the log has a ^in, with — in retarded propagators G, 11 and -|- in advanced 
propagators G*, 11*.) Now we want to compute SIl{Q) using 

flip) = ^V(iVj)x{p) (A6) 

where x{p) = P'^p'^xip)- Then 

^ = V V P ^^^^ ^{i^j)^kX (p) j ■ (A7) 

The integration separates into an angular and a radial part. Integrating the x' radial 
term by parts gives 

dU,,{Q) = ^(^^J pdpxip)^ X A,, = (35mlA,, , 

(7 [ f 2 \ f v^v'^q^ \ 

A^^ ^ J [vi5jk + VjSik - AviVjVk + g^-yfe j [v^g^'' - ^^q^_qo ] ' (^8) 

This depends on Q only through q and q^/q; henceforth we rescale q to be a unit vector, 
and q^ — rj = g°/|q|. 
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First let us find Aqq: 



aij f Vjqj + Vjqj - iVjVjV ■ q + 2SijV ■ q/3 _ aij 



(A9) 



(since this is the only possible tensorial structure). To find A, contract the integral with 
qiqj and define x = q • v: 



A^- dx '— ^2i'irf -1) - Cirf -2r})\TL -. 

4j_i x-T] uj-q 

Next consider A^'^: in practice we will only need 

(r (^ij f ,r> {-Vk + qkq-v)(viqj + Vjqi-AviVjV-q + 25ijV-ci/3) 

2 J v-q-gO 



(AlO) 



= {(Tij/2)B [qiSjk + (Ij^ik - '^mj^lk] ■ 
The coefficient is found by contracting with qiSjk'- 

^ Ilr]-I2r]^ (1 -r/2)(l -4r;2) uj+q 

B — :: 1 : m 



(All) 



(A12) 



6 4 u—q 

Finally we need Ai^. In practice we need it only contracted against transverse projectors: 



7m = — j dfl^ \^ViSjk + VjSik - 4:ViVjVk + '^SijVk^ {^Wmk - ^^'^"^^^0 



A 



(A13) 



must be of form (defining 5im = 5im - qiqm) 



^IrArs^ms 



ClSijn q(iqj) + C2 



25. 



(A14) 



Contracting both the quantity in parenthesis and the original integral expression with two 
independent tensors, such as qtqjSim and SuSjm, determines the coefficients: 



^ (l-772)(157;2-4) 77(1 -772)(3- 5772) co+q 

Gi — \ in 

D 4 LJ—q 

^ (l-^2)(2-3r?2) r]{l-r]y ^^u;+q 



6 



uj—q 



(A15) 



2. Fermionic self-energy 

The fermionic self-energy correction is [39] (convention l/(^ — E)) 



s(g) 



27r2 



J 2p J 



dfip [2/, + /, + /, 



P • 7 - 7 
''p-Cl-q' 



Hence the equilibrium value is 



-2 + ?7 In 



u!+q 
uj—q 



-7° In 



oj+q 
uj—q 



(A16) 



(A17) 
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Taking fi from Eq. (A6), the correction term is 



SE = (3 



1^ J pdp{Xg{p) + Xqip))^ ^ 



-Vklk + 7 
v-q-T) 



= I Pdpixgip) + XM)] A = 



27r2g 

We have A = ^o7° + ^feTfc- 
Start with A°: 



Similarly 



Ak = 



mo) 



X 



1/3 



dx 

I X — 7] 



— 3?7 31]'^ — 1 u+q 



+ 



In 



u—q 



ViVj 



v-q-rj 



a. 



Qi^jk + Qj^ik — -^SijQk 



Determine the coefficients by contracting with qiqjqk and with qiSjk'- 



2ki + Ak2 1 
3 " 6 



dx 



X — 3x^ 

X — Tj 



Therefore 



«1 



K2 = 



4 - Ibr] 
6 

-2 + 377^ 
6 



2ki + 10k2 _ 1 
3 ~ 3 



3// + ^rrf , uj+q 
■ m 



—X 



dx- 

1 X-Tj 



+ - ^In 



uj—q 
oj^-q 
uj—q 



(A18) 



(A19) 



(A20) 



(A21) 



(A22) 



Replacing 7^ — » in Eq. (A16) gives an angular averaged integral and so Q^STi^ = 0, or 
7/A° — (ki + 2/t2) = 0, which is satisfied. This is a fast way to see that the correction to the 
hard propagation velocity is isotropic. 



3. Bosonic 2 <-> 2 contribution to Ci._\4^ 

We work in the plasma rest frame and systematically approximate that the incoming 
particle energies p, k are much larger than the transfer momentum q or frequency < q. 
Using the integration variable parametrization of [19], the contribution to njj^2 order is 

n,„2order 3 J^dp J^dk J^qdq ^dfj J^^^ ^ fo{p)[l±fo{p)] fo{m±Mk)] 

xTxijipf-f [xrsiv) + Xr.(k) - Xr.(p') - X..(k')) [mIM, + M^M\) .(A23) 
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Here Aio,-Mi are to be normalized as in Eq. (A2); we have absorbed all color factors into 
Aah which in SU(A''c) gauge theory with rif fermions is 16dfnfCf/dA for fermion-fermion 
scattering, WdfUfCfCA for fermion-boson scattering and MaCI^ for boson-boson scattering. 
Symmetry between p,p',k,k' allows us to replace 



(A24) 



Xij{p) ^ (xu(p) + Xu(k) - Xij(p') - Xy (k')) 
and small q approximations allow [19] 

(xi,(p)+X»j(k)-Xij(p')-Xij(k')) - -q(3''(2p^iqj)px{p)+VP{iPj)P^x{p)-{p^k)^ (A25) 

and similarly for the xim term. 

All angles are determined by the rj, variables; in particular Xpq = rj = x^q and Xpk = 
rf + {1—rf) COS0. Therefore, extracting a factor oi pk from M.q and M.i, M.q = J^o/pk, 
the integrals over the magnitudes p, k factorize from the integrals over q, u, (j). Defining the 
integrals 



K2 
K, ) 



POD 

= / dpp'^{-M) X < 
Jo 



1 

Ap^XX' 
p'ixT 

2px 

2 — / 

P X 



(A26) 



we need, for the double fermion term for instance. 



-2Klp{iqj)k(^iqm) - 4X4X577P(iPj>k(;q^) - 2i^|77^p(jPj)k(/k^>^ , (A27) 
where we used p k symmetry to simplify some terms. The matrix element squared is 

(A28) 

X ( GlQAqi^rQs) + GooGtB P{rCls) + ^(As) - ^Tjq^As) 



M*oMi = 16{Glo + {l-v'')cos(f)G*j) x /35m^^ 



-Grp 



Ci{l-rf) cos0q(^q^) + 2C2(p(^k^) - vii^i^As) - W{As) + ^7^q(rq5); 



The integral |^ can always be done analytically by replacing cos^^'-*^'^'^'^) = (1, 0, |, 0, |). 
Using repeatedly Eq. (2.44) the evaluation of Eq. (A27) is now straightforward, if lengthy. 

One potential pitfall in performing the g, 7] integrals in Eq. (A27) is the possibility of an 
infrared small q divergence. This can come about because Griq, if) behaves, for q < itl-q and 
77 < q^ /m^, like Gt ~ ^/(f- The integration region over which this behavior applies is q^dq 
but the G*rpG'^ term in Eq. (A28) is 1/q^ so there is a potential log divergence. To determine 
whether this divergence occurs it is sufficient to approximate 77 = in the integrands, other 
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than in Gt- In this hmit Ci = — 2C2. Only the KqKi and Kl terms are zero-order in 77 so 
only they need be computed; the relevant global angular averages are 

KqKi{. . .) COS0 j dQgiobalP(iqi>P(/qm> {v{rK) " COS^q^^q^)) , 

Kl{. . .) COS (f) j rfQgiobaiP(jqj)k(;qr„) (^P(rks) - COS 0q(rqs>) ■ (A29) 

Applying Eq. (2.44) setting Xpq = = x^q and Xpk = coscf) and averaging over 0, one finds 
that each term happens to vanish, so the potential IR divergence does not occur. 



4. Fermionic 2 2 contribution to Ci-Mi 

The infrared region of virtual fermion exchange is also important at leading order for 
transport [19]. The contribution is still described by Eq. (A23) but with A = 32nfCfdf 
each for pair annihilation and Compton scattering. Since the matrix element is less infrared 
singular, we can approximate Xrs{p) — Xrs{p') ^-nd similarly for k. But if x{p) represents 
a fermion, then x{k),x{k') represent a quark and gluon for annihilation, but a gluon and 
quark for Compton scattering. Therefore, summing over the processes, the p, k cross-terms 
cancel and wc may approximate 

Xij{p) (Xrsip) + XrsO^) - Xrs{p') - Xrs{^')) = \ {xijM " XijAP)) {ximM " XlmAP)) 

(A30) 

where the subscripts g, g indicate if the species is a quark or a gluon. This simplifies matters 
considerably; pulling a factor pk out of the p, k integrals we need are 

j dpp'foAp)[^+MPm, - Xof j dkkfoAk)[l+Mk)] (A31) 
which multiply the q, r) integral 

J qdq dii j ^^p^,p^^p^ip^^ [mIM, + h.c.) (A32) 

with 

with = — and ^S^^ as given in Appendix A 2. The trace and global angular 
average are straightforward but tedious. 

5. CoUinear 1 2 contribution to Ci-Mi 

According to [20], the rate at which a particle in the thermal medium splits into two is 
given by 

Ci«2[/(P)] = r dp'dk'5{p-p'-k')^l{p,p',k') 

be 

X (/(p)[l±/(k')][l±/(p')] - [l±/(p)]/(k')/(p')) (A34) 
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where p. k'p' arc collinear at leading order, that is, k' = kp. We saw how this term gives 
rise to contributions to Cu. It also contributes to Ci-Mi because the splitting rate 7^^ is 
sensitive to the details of the plasma, and can be expanded as 

7?c = Co + /i76c,i + ■ ■ ■ (A35) 

We need to evaluate 7^^^; it will then contribute to Ci-Mi through Eq. (A34) with the 
population functions (second line) replaced by 

^ f{k')f{p')[l±f{p)] (3' p'xip) {p'xip) - k"m - p"x{p')) PiiPmiPm) ^ . (A36) 

Besides overall coefficients tabulated in [20], 7^^ is proportional to the integral over the 
solution to an integral equation: 

76» oc J d^h2h-F, 

2h = {iSE)F{h) + J |^C(q^){(a-^)[F(h)-F(h-^V)] 

+ ^[F(h) - F(h -pV)] + ^[F(h) - F(h + pq^)]} . (A37) 

Here h is a vector in the 2-component space transverse to p, 6E is medium dependent but 
in a way which is insensitive to fi (see footnote 8); however C{q±), which represents the 
differential rate to scatter with transverse momentum transfer q^, is sensitive. Explicitly, 

C(qx) = J ^Gl^{q' = q,,q^) (A38) 

with G^_,_ the gauge boson Wightman function, equal to T /uj times the discontinuity in the 
retarded function. In Coulomb gauge this is 

2T 

G%^{Q) = — Disc (Goo + Gt,..) . (A39) 

Here the retarded Green functions Gqo, Gt include the first order corrections, that is, Gt = 
Gt,g-\-Gt,q5TItGt,g- According to [40], analyticity properties allow for the simple evaluation 
of this integral: 

G(qx) = t(GtA^, 0, qx) + Goo(0, 0, q^)) . (A40) 
In equilibrium this reproduces the sum rule of Aurenche, Gelis, and Zaraket [41], 

C(q±) -t{^- -^t^] ■ (A41) 



qi qi + 



For our application the first order shift is 



qi (?i + ^d) 



28 



Using Eq. (A9) , Eq. (AlO), Eq. (A14), and Eq. (A15), 



6Uoo = pSml^q{rqs){-'2) , 



T,zz 



PSml^ (^-'^q/^rqs) + '^P{rPs)^ (A43) 



where we used that the z direction means the p direction. 
U q — X cos + y sin then 

Q{rQs) QrQs ^ '^^rz^sz ~l~ "^{^rx^sx ~l~ ^ry^sy) 

+ - {{drx^sx - SrySsy) COS 20 + {SrxSgy + Sgx^ry) Sm 2(f))^ 

= ~P(rPs} + 0{cos'^(j),sm'^(j)). (A44) 

When expanding Eq. (A37) to hnear order in Ci the (p dependent terms will yield depen- 
dence in F which cancels on angular h integration; therefore these terms may be dropped 
and replaced with —p^rPs}/'^- Hence 

^Hoo = pSml^p^rPs) = SUt,zz (A45) 

To evaluate the shift induced by the correction we have found to Ci, we should expand 
Eq. (A37) linearly in the correction to C(q): schematically (recycling the inner product 
notation for functions over h with J (Ph as inner product) 

|2h) = {iSE + Co + Ci)\F), 

The tools for solving this integral equation are similar to those used in solving the Boltzmann 
equation. The integral we need is (2h | F) . With explicit formulae for everything, the result 
of the analysis is almost straightforward. 

There is one complication, however. Plugging it all in, 

C.(<U) - (KA*,^) (4 + (^^^) (A47) 

has a 1/q^ singularity at small ^j,. Together with the integration measure (Pq,± and the F 
differencing, which on angular averaging behaves like F(h.) — F(h+aq±) ~ a^g^V^-F(h), the 
rest of the integration behaves like q\dq±, resulting in a log IR divergence. The divergence 
is cut off at large momenta by the Debye scale, where F starts to display more complicated 
behavior. In the infrared the calculation becomes unreliable at exchange momentum ~ 
g^T where the perturbative expansion breaks down. We expect that in a nonabelian gauge 
theory the divergence is cut off at this scale, but we are unable to compute the IR end in 
detail. In order to push forward with the calculation we cut the integral off by replacing 
l/q± with l/{q'j_ + (emo)^)^ in the denominator, which allows to extract the coefficient and 
constant under the log. However, the contribution to Ai arising from coUinear contributions 
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to Ci-Mi numerically very small, and the coefficient of this log is still smaller, never 
exceeding 0.003 for 3-flavor QCD and 0.0003 for pure-glue QCD. Therefore in practice the 
uncertainty from resolving this logarithm is too small to see in Figure 2. 
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